{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c9141385",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import integrate\n",
    "import math\n",
    "import numpy as np, matplotlib.pyplot as plt\n",
    "from sklearn.metrics import mean_absolute_error as MAE\n",
    "import time\n",
    "import tqdm\n",
    "from tqdm import trange\n",
    "from environment import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7520103a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "Liu_all=pd.read_csv('Liu_all.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "66f7e028",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Liu_all['ar_ratio']=np.random.uniform(1.3,1.8,len(Liu_all))\n",
    "# Liu_all.to_csv('Liu_all.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8e358c83",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "env= environment()\n",
    "jc_scale=6\n",
    "env.update_fitted_para(0.01,0.03)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ba919dfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "tcd8_target=Liu_all['tcd8_pop']\n",
    "r_l=np.ones(len(Liu_all))*0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "16d621e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def opt_fix_para():\n",
    "    time_l=[]\n",
    "    cd8_l=[]\n",
    "    for i in trange(len(Liu_all)):\n",
    "        mu=Liu_all.loc[i,'mu']\n",
    "        r=r_l[i]\n",
    "        jc=Liu_all.loc[i,'mhc_norm']\n",
    "        jc_=(1-np.exp(-jc*jc_scale))/(1+np.exp(-jc*jc_scale))\n",
    "        R_L_prop=Liu_all.loc[i,'Treg_pro']\n",
    "        tumor=Liu_all.loc[i,'tumor_size']\n",
    "        a=Liu_all.loc[i,'a']\n",
    "        ar_ratio=Liu_all.loc[i,'ar_ratio']\n",
    "        t1=Liu_all.loc[i,'time']\n",
    "        \n",
    "        x_=env.run(mu,r,jc_,R_L_prop,t1,a,ar_ratio)\n",
    "        j=cal_time(x_,tumor)\n",
    "        cd8_pred=min(x_[j][2]/(x_[j][0]+x_[j][1]),1)\n",
    "        cd8_l.append(cd8_pred)\n",
    "#         print(j/24,x_[j])\n",
    "        time_l.append(j/24+10)\n",
    "    return MAE(cd8_l,tcd8_target),cd8_l,time_l\n",
    "def cal_time(x_,tumor):\n",
    "    for i in range(len(x_.T[0])):\n",
    "        if x_.T[0][i]+x_.T[1][i]>tumor:\n",
    "            break\n",
    "    return i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a02d2717",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████| 120/120 [02:03<00:00,  1.03s/it]\n"
     ]
    }
   ],
   "source": [
    "#initial r=0.1\n",
    "Liu_all['time']=int(600)\n",
    "mae,cd8_l,time_l=opt_fix_para()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "20c2105b",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.011121013692878202"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mae"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "5f9d6857",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_l_=[]\n",
    "for i in time_l:\n",
    "    time_l_.append(int(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "a0d20b89",
   "metadata": {},
   "outputs": [],
   "source": [
    "Liu_all['time']=time_l_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "c699805f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████| 120/120 [05:37<00:00,  2.81s/it]\n"
     ]
    }
   ],
   "source": [
    "#fit r\n",
    "for i in trange(len(Liu_all)):\n",
    "    tcd8_target_s=tcd8_target[i]\n",
    "    r_origin=r_l[i]\n",
    "    ar_ratio=Liu_all['ar_ratio'][i]\n",
    "    error_origin=np.abs(tcd8_target_s-cd8_l[i])\n",
    "    if tcd8_target_s>cd8_l[i]:\n",
    "        r=up_adjust(i,tcd8_target_s,r_origin,error_origin,ar_ratio)\n",
    "    else:\n",
    "        r=down_adjust(i,tcd8_target_s,r_origin,error_origin,ar_ratio)\n",
    "    r_l[i]=r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "c6fdecc4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████| 120/120 [00:45<00:00,  2.64it/s]\n"
     ]
    }
   ],
   "source": [
    "#calculate mae after r fitting\n",
    "mae,cd8_l,time_l=opt_fix_para()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "b362c9a1",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.001484007401008113"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#mae after r fitting\n",
    "mae"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "e6c25ab4",
   "metadata": {},
   "outputs": [],
   "source": [
    "jc_scale=6\n",
    "\n",
    "def down_adjust(i,cd8_target,r_origin,error_origin,ar_ratio):\n",
    "    mu=Liu_all.loc[i,'mu']\n",
    "    r_last=r_origin\n",
    "    error_last=error_origin\n",
    "    error=error_origin\n",
    "    while error<=error_last:\n",
    "        r=r_last-0.002\n",
    "        jc=Liu_all.loc[i,'mhc_norm']\n",
    "        jc_=(1-np.exp(-jc*jc_scale))/(1+np.exp(-jc*jc_scale))\n",
    "        R_L_prop=Liu_all.loc[i,'Treg_pro']\n",
    "        tumor=Liu_all.loc[i,'tumor_size']\n",
    "        a=Liu_all.loc[i,'a']\n",
    "\n",
    "        t1=Liu_all.loc[i,'time']\n",
    "        x_=env.run(mu,r,jc_,R_L_prop,t1,a,ar_ratio)\n",
    "        j=cal_time(x_,tumor)\n",
    "    #     print('cd8+T:',x_[i][2])\n",
    "        cd8_pred=min(x_[j][2]/(x_[j][0]+x_[j][1]),1)\n",
    "        error_last=error\n",
    "        error=np.abs(cd8_pred-cd8_target)\n",
    "        r_last=r\n",
    "    return r+0.002\n",
    "\n",
    "def up_adjust(i,cd8_target,r_origin,error_origin,ar_ratio):\n",
    "    mu=Liu_all.loc[i,'mu']\n",
    "    r_last=r_origin\n",
    "    error_last=error_origin\n",
    "    error=error_origin\n",
    "    while error<=error_last:\n",
    "        r=r_last+0.002\n",
    "        jc=Liu_all.loc[i,'mhc_norm']\n",
    "        jc_=(1-np.exp(-jc*jc_scale))/(1+np.exp(-jc*jc_scale))\n",
    "        R_L_prop=Liu_all.loc[i,'Treg_pro']\n",
    "        tumor=Liu_all.loc[i,'tumor_size']\n",
    "        a=Liu_all.loc[i,'a']\n",
    "        \n",
    "        t1=Liu_all.loc[i,'time']\n",
    "        x_=env.run(mu,r,jc_,R_L_prop,t1,a,ar_ratio)\n",
    "        j=cal_time(x_,tumor)\n",
    "    #     print('cd8+T:',x_[i][2])\n",
    "        cd8_pred=min(x_[j][2]/(x_[j][0]+x_[j][1]),1)\n",
    "        error_last=error\n",
    "        error=np.abs(cd8_pred-cd8_target)\n",
    "        r_last=r\n",
    "        if r>0.25:\n",
    "            break\n",
    "    return r-0.002   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "0b2cf264",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate error for jl fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "1cdf9d71",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def cal_jl_error(jl):\n",
    "\n",
    "    cycle_l=[]\n",
    "    count=0\n",
    "    target=[]\n",
    "    pred=[]\n",
    "    for i in trange(120):\n",
    "        target_re=Liu_all.loc[i,'BR']\n",
    "        if  target_re=='MR' or target_re=='SD':\n",
    "            pass\n",
    "        else:\n",
    "            mu=Liu_all.loc[i,'mu']\n",
    "            r=r_l[i]\n",
    "\n",
    "            jc=Liu_all.loc[i,'mhc_norm']\n",
    "            jc_=(1-np.exp(-jc*jc_scale))/(1+np.exp(-jc*jc_scale))\n",
    "            R_L_prop=Liu_all.loc[i,'Treg_pro']\n",
    "            tumor=Liu_all.loc[i,'tumor_size']\n",
    "            a=Liu_all.loc[i,'a']\n",
    "            ar_ratio=Liu_all.loc[i,'ar_ratio']\n",
    "\n",
    "            t1=Liu_all.loc[i,'time']\n",
    "            x_=env.run(mu,r,jc_,R_L_prop,t1,a,ar_ratio)\n",
    "            j=cal_time(x_,tumor)\n",
    "            tumor_size_start=x_[j][0]+x_[j][1]\n",
    "            cycles=Liu_all.loc[i,'cyclesOnTherapy']\n",
    "            try:\n",
    "                cycles=int(cycles)\n",
    "            except:\n",
    "                cycles=6\n",
    "\n",
    "            x_begin=x_[j].tolist()\n",
    "            x_begin.append(1)\n",
    "            cd8_prop=Liu_all.loc[i,'tcd8_pop']\n",
    "    #         print(cd8_prop,R_L_prop)\n",
    "            x_begin[2]=cd8_prop*tumor_size_start\n",
    "            for c in range(cycles):\n",
    "                x_begin[-1]=1\n",
    "                x_end=env.step(x_begin,mu,r,jc_,R_L_prop,jl)\n",
    "                x_begin=x_end\n",
    "    #             if x_end[0]+x_end[1]<1e2:\n",
    "    #                 cycle_l.append(c)\n",
    "    #                 break\n",
    "    #             print(x_begin)\n",
    "            tumor_size_end=x_end[0]+x_end[1]\n",
    "            pred_re=1-tumor_size_end/tumor_size_start\n",
    "            if pred_re>0:\n",
    "                \n",
    "                pred.append(1)\n",
    "    #               \n",
    "            else:\n",
    "                pred.append(0)\n",
    "\n",
    "\n",
    "            if target_re=='PR' or target_re=='CR':\n",
    "                target.append(1)\n",
    "            else:\n",
    "                target.append(0)\n",
    "\n",
    "\n",
    "    \n",
    "    return target,pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "a755179b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score,recall_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "a4fc5916",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████| 120/120 [01:21<00:00,  1.47it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6310679611650486\n",
      "0.6382978723404256\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "target,pred_2=cal_jl_error(0.11)\n",
    "print(accuracy_score(target,pred_2))\n",
    "print(recall_score(target,pred_2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "a950c455",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████| 120/120 [01:20<00:00,  1.48it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6407766990291263\n",
      "0.574468085106383\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "target,pred_1=cal_jl_error(0.1)\n",
    "print(accuracy_score(target,pred_1))\n",
    "print(recall_score(target,pred_1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "787a7a2f",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████| 120/120 [01:21<00:00,  1.48it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6213592233009708\n",
      "0.723404255319149\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "target,pred_3=cal_jl_error(0.12)\n",
    "print(accuracy_score(target,pred_3))\n",
    "print(recall_score(target,pred_3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "714cb223",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: 64, 1: 13, -1: 26}"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count={0:0,1:0,-1:0}\n",
    "for i in np.array(target)-np.array(pred_3):\n",
    "    count[i]+=1\n",
    "count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "d599d789",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: 66, 1: 20, -1: 17}"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count={0:0,1:0,-1:0}\n",
    "for i in np.array(target)-np.array(pred_1):\n",
    "    count[i]+=1\n",
    "count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "cfb4cbe5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: 65, 1: 17, -1: 21}"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count={0:0,1:0,-1:0}\n",
    "for i in np.array(target)-np.array(pred_2):\n",
    "    count[i]+=1\n",
    "count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f58d1cf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
